
clearvars
starttime = datetime('now');

% ////////////////////////////////////////////////////////
% SET UP FILES AND PARAMETERS
% ////////////////////////////////////////////////////////

main_base = " "; % designate a base folder

bname = " "; % name this set of simulations


nrun = 10;% number of runs 

% set up file names to store data 
fnames = strings(1,nrun);
for i=1:nrun
    fnames(i) = strcat(main_base, "/",bname, "_", string(i));
end

% specify parameters

% [1VISC, 2PACKING, 3BORD/CORN removal, 4REP STR, 5ATT STR, 6ANGLE, 7ACTIVE WIDTH,
% 8ACTIVE SPEED 9CAPT RECT 10ATTR ARROWS 11ML ATTR MULTIPILER, 12ML CAPTURE MULTIPLIER, 13 PLAYGROUND SPEED]; %


%contraction sample parameters
ps(:,1:10) = repmat([30, 1.0, 1, 1,0.15, 25, 0, 0, 0,0, 6, 1, 0.003*10000]',1,10); % base of 0827

% crawling sample parameters
ps(:,11:20) = repmat([30, 1.0, 1, 1,0.15, 0, 50, 4, 0,0, 1, 1, 0.003*10000]',1,10); %*BP

% capture sample parameters
ps(:,21:30) = repmat([30, 1.0, 1, 1,0.15, 0, 0, 0, 100,0, 1, 8, 0.003*10000]',1,10); %*BP


% run and analyze all sims
for q=1:nrun %:nrun  %20:30  %1:nrun
    javaaddpath 'C:\Program Files\MATLAB\R2019b\java\mij.jar' % set Matlab-Image J path
    javaaddpath 'C:\Program Files\MATLAB\R2019b\java\ij.jar' % set Image J path
    MIJ.start('C:\Users\Lab\Fiji.app')
    basedir = fnames(q);
    
    % remove duplicate files if re-running with same base file name 
    if exist(basedir, 'dir')
        a = rmdir(basedir, 's');
    end
    
    mkdir(basedir)
    mkdir(strcat(basedir, "/rois"));
    
    % write base file for sim to read
    bstore = strcat('C:\Users\Lab\Downloads\basefile.txt');
    fid=fopen(bstore,'wt');
    fprintf(fid, '%s', basedir');
    fclose(fid);
    
    % write parameter file for sim to read
    pstore = strcat(basedir, '\params.txt');
    fid=fopen(pstore,'wt');
    fprintf(fid, '%f \n', ps(:,q)');
    fclose(fid);
    
    % current set of parameters
    ps_now = ps(:,q);
    
    
    % ////////////////////////////////////////////////////////
    % MACRO RUNNING
    % ////////////////////////////////////////////////////////
    
    cd(" "); % folder with strain codes 
    
    MIJ.run("Install...", "install=[ ]") % file path for main sim code (CE_sim_SA.ijm)
    MIJ.run("Run CE");
    
    %COMMENT THIS OUT IF RUNNING ALL ANALYSES
%     MIJ.run("Clear Results");
%     MIJ.closeAllWindows
%     MIJ.exit
%     javarmpath 'C:\Program Files\MATLAB\R2019b\java\mij.jar'
%     javarmpath 'C:\Program Files\MATLAB\R2019b\java\ij.jar'
    
    %%{
    
   % //comment block start for just running sim and segmentation, no analysis
    %these need to be installed to the local FIJI as plugins
    MIJ.run("cellstrain new box");

% ////////////////////////////////////////////////////////
% ANALYSIS
% ////////////////////////////////////////////////////////
    
% ********* cell and corona strain analysis **************
    
    % Import Cellstrain

    opts = delimitedTextImportOptions("NumVariables", 4);
    
    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = "\t";
    
    % Specify column names and types
    opts.VariableNames = ["SliceNumA", "CellID", "CentroidX", "CentroidY", "Angle", "major", "minor", "Area", "Width", "Height", "Perim"];
    opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
    
    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";
    
    % Import the data
    CellEllipse = readtable(strcat(basedir, "/CellEllipse.txt"), opts);
    
    
    % Clear temporary variables
    clear opts
    
     if (isa(CellEllipse,'table'))
        CellEllipse = table2struct(CellEllipse);
     end
    
     
     times =  [CellEllipse.SliceNumA];
     time = times(length(times));
    
    tr=[CellEllipse([CellEllipse.SliceNumA]==floor(time/2)).CellID];
    tr=unique(tr);
    a=[CellEllipse([CellEllipse.SliceNumA]==1).CellID];
    b=[CellEllipse([CellEllipse.SliceNumA]==time).CellID];
     c=setdiff(b,a);
     [bn,ibn,ib]=intersect(tr,c);
     tr(ibn)=[];
     
    fid=fopen('C:\Users\Lab\Downloads\tr.txt','wt');
    fprintf(fid, '%.0f \n', tr);
    fclose(fid);
     
     
    MIJ.run("corona one level"); 
    MIJ.run("t1track");
    
    MIJ.run("Clear Results");
    MIJ.closeAllWindows
    MIJ.exit
    javarmpath 'C:\Program Files\MATLAB\R2019b\java\mij.jar'
    javarmpath 'C:\Program Files\MATLAB\R2019b\java\ij.jar'
    
    % Import CoronaData

    % Setup the Import Options and import the data
    opts = delimitedTextImportOptions("NumVariables", 12);
    
    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = "\t";
    
    % Specify column names and types
    opts.VariableNames = ["CellID", "SliceNumA", "CoronaCells", "CentroidX", "CentroidY", "Angle", "major", "minor", "Width", "Height", "Area", "sided"];
    opts.VariableTypes = ["double", "double", "string", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
    
    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";
    
    % Specify variable properties
    opts = setvaropts(opts, "CoronaCells", "WhitespaceRule", "preserve");
    opts = setvaropts(opts, "CoronaCells", "EmptyFieldRule", "auto");
    
    % Import the data
    CoronaData = readtable(strcat(basedir, "\Coronas\CoronaData.txt"), opts);
    
    
    % Clear temporary variables
    clear opts
    
    % Setup the Import Options and import the data
    opts = delimitedTextImportOptions("NumVariables", 1);
    
    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = ",";
    
    % Specify column names and types
    opts.VariableNames = "Value";
    opts.VariableTypes = "double";
    
    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";
    
    
    % Clear temporary variables
    clear opts
    %
    cd(" "); % return to folder with strain analysis code 
    
    
    
    cells = unique([CoronaData.CellID]);
    
    % endpoint strain calculation
    [CellStrain_end,~] = CellStrainCalc_endpoint(CellEllipse, time, cells);
    [CoronaStrain_end,~] = CellStrainCalc_endpoint(CoronaData,time, cells);
    
    % step strain calculation
    [CellStrain,~] = CellStrainCalc(CellEllipse,cells, time);
    [CoronaStrain,~] = CellStrainCalc(CoronaData,cells,time);
    
   
   
    % Get data out
    if (isa(CoronaData,'table'))
        CoronaData = table2struct(CoronaData);
    end

   
    cd(" "); % return to folder with strain analysis code 

    
    G = zeros(time,length(cells));
    GG = zeros(time,length(cells));
    K = zeros(time,length(cells));
    KK = zeros(time,length(cells));
    M = zeros(time,length(cells));
    MM = zeros(time,length(cells));
    g = zeros(time, 1);
    gg = zeros(time, 1);
    k = zeros(time, 1);
    kk = zeros(time, 1);
    m = zeros(time, 1);
    mm = zeros(time, 1);
    
    
    for jj=1:length(cells)
        % g is corona strain in the xx direction.
        % gg is cell strain in the xx direction
        
        % k is corona strain Eyy
        % kk is cell strain in the yy direction
        % cells(jj)
        
        % m is Area strain of corona
        % mm is Area strain of cell
        
        for i=2:time
            
            %collect Exx of that cell's corona
            g(i)=[CoronaStrain_end([CoronaStrain_end.CellID]==cells(jj) & [CoronaStrain_end.time]==i).Exx];
            
            
            %Doing the same as above for the YY direction:
            k(i)=[CoronaStrain_end([CoronaStrain_end.CellID]==cells(jj) & [CoronaStrain_end.time]==i).Eyy];
            
            
            gg(i)=[CellStrain_end([CellStrain_end.CellID]==cells(jj)& [CellStrain_end.time]==i).Exx];
            
            kk(i)=[CellStrain_end([CellStrain_end.CellID]==cells(jj)& [CellStrain_end.time]==i).Eyy];
            
            %AC is area strain
            mm(i)= [CellStrain_end([CellStrain_end.CellID]==cells(jj)& [CellStrain_end.time]==i).EA];
            m(i) = [CoronaStrain_end([CoronaStrain_end.CellID]==cells(jj)& [CoronaStrain_end.time]==i).EA];
            
        end
        
        G(:,jj) = g;
        GG(:,jj) = gg;
        K(:,jj) = k;
        KK(:,jj) = kk;
        M(:,jj) = m;
        MM(:,jj) = mm;
        
    end
    
    S = zeros(time,length(cells));
    SS = zeros(time,length(cells));
    V = zeros(time,length(cells));
    VV = zeros(time,length(cells));
    U = zeros(time,length(cells));
    UU = zeros(time,length(cells));
    s = zeros(time, 1);
    ss = zeros(time, 1);
    v = zeros(time, 1);
    vv = zeros(time, 1);
    u = zeros(time, 1);
    uu = zeros(time, 1);
    
    for jj=1:length(cells)
        % g is corona strain in the xx direction.
        % gg is cell strain in the xx direction
        
        % k is corona strain Eyy
        % kk is cell strain in the yy direction
        %cells(jj)
        
        % m is Area strain of corona
        % mm is Area strain of cell
        
        for i=2:time
            %   i
            %collect Exx of that cell's corona
            s(i)=[CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Exx];
            
            
            %Doing the same as above for the YY direction:
            v(i)=[CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Eyy];
            
            
            ss(i)=[CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Exx];
            
            vv(i)=[CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Eyy];
            
            %AC is area strain stuff
            uu(i)= [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).EA];
            u(i) = [CoronaStrain([CoronaStrain.CellID]==cells(jj)& [CoronaStrain.time]==i).EA];
            
        end
        
        S(:,jj) = s;
        SS(:,jj) = ss;
        V(:,jj) = v;
        VV(:,jj) = vv;
        U(:,jj) = u;
        UU(:,jj) = uu;
        
    end
    

% ********* Markov analysis **************

    concordx = zeros(length(cells),time);
    concordy = zeros(length(cells),time);
 
    
    % Finding concordance for Exx:
    % Iterate over cells and time to populate the concordx matrix
    for jj = 1:length(cells)
        
        for i = 1:time
            %  Find Exx for the Corona for the cell that the loop is on, for
            %  the timestep that the loop is on
            Cor = [CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Exx];
            
            % Find Exx for the cell & time that the loop is on
            Cell = [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Exx];
            
            % Check if the sign of Exx for the corona is the same as the cell
            % in question at time t
            if sign(Cor) == sign(Cell)
                % If the sign is the same, call this concordance and assign
                % value 1
                concordx(jj,i) = 1;
            else
                % Let 0 signify discordance for that cell and time if the sign
                % does not match
                concordx(jj,i) = 0;
            end
        end
        
    end
    
    
    % find the concordance for Eyy with new matrix name
    
 
    % Iterate over cells and time to populate the concord1 matrix
    for jj = 1:length(cells)
        
        for i = 1:time
            %  Find Exx for the Corona for the cell that the loop is on, for
            %  the timestep that the loop is on
            Cor = [CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Eyy];
            
            % Find Exx for the cell & time that the loop is on
            Cell = [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Eyy];
            
            % Check if the sign of Exx for the corona is the same as the cell
            % in question at time t
            if sign(Cor) == sign(Cell)
                % If the sign is the same, call this concordance and assign
                % value 1
                concordy(jj,i) = 1;
            else
                % Let 0 signify discordance for that cell and time if the sign
                % does not match
                concordy(jj,i) = 0;
            end
        end
        
    end
    
    P1x =  zeros(length(cells),time);
    P2x =  zeros(length(cells),time);
    P3x =  zeros(length(cells),time);
    P4x =  zeros(length(cells),time);
    
    for j = 1:length(cells)
        y_obs = concordx(j,:)';
        y_obs(y_obs == 1) = 2;
        y_obs(y_obs < 1) = 1;
        L = length(y_obs);
        T = L;
        
        P_MC = zeros(2,2);
        
        for t=1:L-1
            P_MC(y_obs(t),y_obs(t+1))= P_MC(y_obs(t),y_obs(t+1))+1;
            
            P1x(j,t) =  P_MC(1,1)/sum(P_MC(1,:));
            P2x(j,t) =  P_MC(1,2)/sum(P_MC(1,:));
            P3x(j,t) =  P_MC(2,1)/sum(P_MC(2,:));
            P4x(j,t) =  P_MC(2,2)/sum(P_MC(2,:));
        end
        
    end
    
    P1y =  zeros(length(cells),time);
    P2y =  zeros(length(cells),time);
    P3y =  zeros(length(cells),time);
    P4y =  zeros(length(cells),time);
    
    for j = 1:length(cells)
        y_obs = concordy(j,:)';
        y_obs(y_obs == 1) = 2;
        y_obs(y_obs < 1) = 1;
        L = length(y_obs);
        T = L;
        
        P_MC = zeros(2,2);
        
        for t=1:L-1
            P_MC(y_obs(t),y_obs(t+1))= P_MC(y_obs(t),y_obs(t+1))+1;
            
            P1y(j,t) =  P_MC(1,1)/sum(P_MC(1,:));
            P2y(j,t) =  P_MC(1,2)/sum(P_MC(1,:));
            P3y(j,t) =  P_MC(2,1)/sum(P_MC(2,:));
            P4y(j,t) =  P_MC(2,2)/sum(P_MC(2,:));
        end
        
    end
    
% ********* signed strain analysis **************
    
    sstrainx = zeros(length(cells),time);

    for jj = 1:length(cells)
        
        for i = 2:time
            %  Find Exx for the Corona for the cell that the loop is on, for
            %  the timestep that the loop is on
            Cor = [CoronaStrain([CoronaStrain_end.CellID]==cells(jj) & [CoronaStrain.time]==i).Exx];
            
            % Find Exx for the cell & time that the loop is on
            Cell = [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Exx];
            
            % Check if the sign of Exx for the corona is the same as the cell
            % in question at time t
            if sign(Cor) == sign(Cell)
                % If the sign is the same, call this concordance and assign
                % value 1
                sstrainx(jj,i) = abs(Cell - Cor);
                
            else
                % Let 0 signify discordance for that cell and time if the sign
                % does not match
                sstrainx(jj,i) = -(abs(Cell-Cor)) ;
            end
        end
        
    end
   
    
    sstrainy = zeros(length(cells),time);

    for jj = 1:length(cells)
        
        for i = 2:time
            %  Find Exx for the Corona for the cell that the loop is on, for
            %  the timestep that the loop is on
            Cor = [CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Exx];
            
            % Find Exx for the cell & time that the loop is on
            Cell = [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Exx];
            
            % Check if the sign of Exx for the corona is the same as the cell
            % in question at time t
            if sign(Cor) == sign(Cell)
                % If the sign is the same, call this concordance and assign
                % value 1
                sstrainy(jj,i) = abs(Cell - Cor);
                
            else
                % Let 0 signify discordance for that cell and time if the sign
                % does not match
                sstrainy(jj,i) = - ( abs(Cell-Cor) );
            end
        end
        
    end
    
    
% ********* Shape index and area analyis **************
    
    %CellEllipse = table2struct(CellEllipse);
    
    si = zeros(length(cells), time);
    areas = zeros(length(cells), time);
    ars = zeros(length(cells), time);
    angles = zeros(length(cells), time);
    
    for j=1:length(cells)
        for i=1:time
            perim = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).Perim;
            area = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).Area;
            si(j,i) = perim/sqrt(area);
        end
    end
    
     for j=1:length(cells)
        for i=1:time
            
            area = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).Area;
            areas(j,i) = area;
            ang = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).Angle;
            angles(j,i) = ang;
            
            major = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).major;
            minor = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).minor;
            
            ars(j,i) = major/minor;
            
        end
     end
    

         

% ********* MSD analysis **************

    % Use x and y centroids from CellEllipse
    
    
    % Calculate MSDs
    
    % remove first zero, create a new matrix with each cell and time from
    % long array
    
    xCens = zeros(length(cells), time);
    yCens = zeros(length(cells), time);
    
    for i=1:length(cells)
        for j=1:time
            xCens(i,j)=CellEllipse([CellEllipse.SliceNumA]==j & [CellEllipse.CellID]==cells(i)).CentroidX;
            yCens(i,j)=CellEllipse([CellEllipse.SliceNumA]==j & [CellEllipse.CellID]==cells(i)).CentroidY;
        end
    end
    
    
    MSD = zeros(length(cells), time);
    
    for i = 1:length(cells) % over cells
        for j = 1:time  % over time
            sums(i,j) =  ( ( xCens(i,j) - xCens(i,1) )^2 + (  yCens(i,j) -  yCens(i,1) )^2 );
        end
    end
    
      cs = cumsum(sums,2);
    asum = cumsum(areas,2); 

    for i = 1:length(cells)
        for j = 1:time
            MSD(i,j) = ( (1/j) .* cs(i,j) )/(asum(i,j) );
        end
    end
    
    
% ********* T1 analysis **************
    
    % import data
    % Setup the Import Options and import the data
    opts = delimitedTextImportOptions("NumVariables", 6);

    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = ["\t", ","];

    % Specify column names and types
    opts.VariableNames = ["t1time", "t1cells", "VarName3", "VarName4", "VarName5", "VarName6"];
    opts.VariableTypes = ["double", "double", "double", "double", "double", "double"];

    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";

    % Import the data
    t1s = readtable(strcat(basedir,"\t1s.txt"), opts);

    % Convert to output type
    t1s = table2array(t1s);

    % Clear temporary variables
    clear opts
    
    %remove duplicate entries, get rid of NaNs
    t1s(isnan(t1s))=0;
    t1s = unique(t1s,'rows');
    
    %create table of binary events by cell and time
    t1 = zeros(length(cells), time);
    
    %t1 = zeros(length(t1s), time); for testing
    for i=1:size(t1s,1)
        for j=2:size(t1s,2)
            if ismember(t1s(i,j), cells)
                cell_ind = find(cells==t1s(i,j));
                %cell_ind = t1s(i,j); for testing
                t1(cell_ind, t1s(i,1)) = 1;
                
            end
        end
    end
    
% ********* SIDEDNESS ANALYSIS **************
    
    %doing initial sidedness of cell by finding n cells in one level corona
    %then doing for other frames, added onto corona macro
    
    sided = zeros(length(cells), time);
    %CellEllipse([CellEllipse.SliceNumA]==j & [CellEllipse.CellID]==cells(i)).CentroidX;
    for j=1:length(cells)
        for i=1:time
            sided(j,i) = CoronaData([CoronaData.CellID]==cells(j) & [CoronaData.SliceNumA]==i).sided;
        end
    end
    
        
% ********* AP STRAIN **************

    % Setup the Import Options and import the data
    opts = delimitedTextImportOptions("NumVariables", 1);

    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = ",";

    % Specify column names and types
    opts.VariableNames = "Value";
    opts.VariableTypes = "double";

    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";

    % Import the data
    apstrain = readtable(strcat(basedir,"\apstrain.txt"), opts);

    % Convert to output type
    apstrain = table2array(apstrain);

% Clear temporary variables
    clear opts
    %load file and populate it for all cells in that field
        apstrains = zeros(length(cells), time);
    %CellEllipse([CellEllipse.SliceNumA]==j & [CellEllipse.CellID]==cells(i)).CentroidX;
    for j=1:length(cells)
        for i=1:time
            apstrains(j,i) = apstrain(i);
        end
    end
    
% stopper/time to 15%

opts = delimitedTextImportOptions("NumVariables", 1);

% Specify range and delimiter
opts.DataLines = [1, Inf];
opts.Delimiter = ",";

% Specify column names and types
opts.VariableNames = "VarName1";
opts.VariableTypes = "double";

% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% Import the data
stopper = readtable(strcat(basedir,"\stopper.txt"), opts);

% Convert to output type
stopper = table2array(stopper);

% Clear temporary variables
clear opts

stopts = zeros(length(cells), time);

        for j=1:length(cells)
            for i=1:time
            stopts(j,i) = stopper;
            end
        end

    
% ////////////////////////////////////////////////////////
% GENERATE 3D MATRIX WITH VARIABLES OF INTEREST
% ////////////////////////////////////////////////////////
    
    % first index is data, second is cell, third is time
    
    data = [];
        
    for i=1:length(cells)
        for j=1:time
            
            data(1,i,j) = cells(i); % cell ID
            
            % populate cell and corona strains
            data(2,i,j) = GG(j,i); % cell exx
            data(3,i,j) = KK(j,i); % cell eyy
            data(4,i,j) = G(j,i); % cor exx
            data(5,i,j) = K(j,i); % cor eyy
            data(6,i,j) = MM(j,i); % cell e_area
            data(7,i,j) = M(j,i); % cor e_area
            
            % populate state transition probabilities
            data(8,i,j) = P1x(i,j); % pdd x
            data(9,i,j) = P2x(i,j); % pdc x
            data(10,i,j) = P3x(i,j); % pcd x
            data(11,i,j) = P4x(i,j); % pcc x
            data(12,i,j) = P1y(i,j); % pdd y
            data(13,i,j) = P2y(i,j); % pdc y
            data(14,i,j) = P3y(i,j); % pcd y
            data(15,i,j) = P4y(i,j); % pcc y
            
            % populate signed difference in cell and corona strain
            data(16,i,j) = sstrainx(i,j); % cell-cor x
            data(17,i,j) = sstrainy(i,j); % cell-cor y
           % add area after com
            
            %populate shape index
            data(18,i,j) = si(i,j);
            
            %populate MSD
            data(19,i,j) = MSD(i,j);
            
            %populate T1s
            data(20,i,j) = t1(i,j);
            
            %populate parameters
            data(21,i,j) = ps_now(1); %viscosity
            data(22,i,j) = ps_now(2); %packing
            data(22,i,j) = ps_now(3); %border/corner
            data(23,i,j) = ps_now(4); %repulsion strength
            data(24,i,j) = ps_now(5); %attraction strength
            data(25,i,j) = ps_now(6); %aniso attrac (angle range)
            data(26,i,j) = ps_now(7); %width of active midline
            data(27,i,j) = ps_now(8); %speed of active midline
            data(28,i,j) = ps_now(9); %attract in rect--> capture width
           
           
           %STEP STRAINS
            data(29,i,j) = SS(j,i); % cell exx
            data(30,i,j) = VV(j,i); % cell eyy
            data(31,i,j) = S(j,i); % cor exx
            data(32,i,j) = V(j,i); % cor eyy
            data(33,i,j) = UU(j,i); % cell e_area
            data(34,i,j) = U(j,i); % cor e_area
            
            
                data(35,i,j) = xCens(i,j);% xc
            data(36, i, j) = yCens(i,j);% yc
            data(37,i,j) = areas(i,j); % area cell
           data(38, i ,j) = ars(i,j);% aspect ratio major/minor
           data(39,i,j)= angles(i,j);% angle
           data(40,i,j) = sided(i,j); %sidedness
            
            %more parameters
            data(41,i,j) = ps_now(10); %attr arrows
            data(42,i,j) = ps_now(11); %attr bias multiplier
            data(43,i,j) = ps_now(12); % capture ML multiplier
            data(44,i,j) = apstrains(i,j); % ap strain from stopper
            data(45,i,j) = ps_now(13); % playground rate
            data(46,i,j) = stopts(i,j); % ts that stopper kicked in 
            
           
        end
    end
    
    cd(basedir);
    save('wspace.mat');
    
    %}
    
    %comment block end for just running sim and segmentation, no analysis
    
    
    % ////////////////////////////////////////////////////////
    % POST ANALYSIS END MATTER
    % ////////////////////////////////////////////////////////
    
    
    
    clearvars -except fnames starttime ps notes
    
end
endtime = datetime('now');
endtime-starttime